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In this paper we report the results of simulations of a 2D gravity driven, dissipative granular 
flow through a hopper system. Measurements of impulse distributions P(I) on the simulated sys- 
tem show flow-velocity-invariant behavior of the distribution for impulses larger than the average 
impulse < I >. For small impulses, however, P(I) decreases significantly with flow velocity, a 
phenomenon which can be attributed exclusively to collisions between grains undergoing frequent 
collisions. Visualizations of the system also show that these frequently colliding particles tend to 
form increasingly large linear clusters as the flow velocity decreases. A model is proposed for the 
form of P(I), given distributions of cluster size and velocity, which accurately predicts the observed 
form of the distribution. Thus the impulse distribution provides some insight into the formation 
and properties of these "dynamic" force chains. 



Introduction Granular materials exhibit a wide spec- 
trum of behavior ranging from gaseous to liquid to solid. 
Remarkably, all of these phases of granular matter re- 
spond to external stimuli in a manner strikingly different 
from ordinary fluids and solids 0, 0|. In static granu- 
lar piles, the spatially inhomogeneous manner in which 
stress is transmitted from the bulk of the pile to the 
boundary has been well documented^ 0]. Experiments 
have shown that the force distribution P(f) at the walls 
is exponential at large forces and exhibits a plateau at 
small forces 5]. In addition, the highly stressed grains 
in static packings are organized into linear structures 
termed "force chains" 0, 0] . The appearance of similar 
large scale structures in flowing granular matter would 
have significant implications for both continuum theories 
and descriptions of jamming in non-thermal systems. A 
recent proposal for a unified picture of jamming in ther- 
mal and non-thermal systems suggests that jamming oc- 
curs due the formation of force chains whose presence is 
signalled by the appearance of a plateau in P(f) @. A 
continuum description of steady-state flow on an inclined 
plane models the system as a collection of transient ID 
chains immersed in a viscous fluidQ. Indeed, transient 
"clusters" have been identified experimentally in granu- 
lar surface flows 8] and shear flows 0, and simulations 
of chute flows have shown evidence for a plateau in P(f) 
as the system approaches iamming|lOj. 

Recent experiments have been performed in a two di- 
mensional hopper geometry to explore the presence of 
incipient force chains in purely collisional gravity-driven 
flow Measurements of the time trace of the im- 

pulse delivered to a transducer placed at the side wall of 
the hopper have shown that the distribution of impulses, 
P(I), displays an exponential decay at large /, as for the 
case of static materials. This exponential form of the 
distribution is maintained for all flow velocities from the 
largest measured (vf = 60.0 cm/s) to the minimum flow 
velocity prior to the point at which the system no longer 
exhibits sustained flow (vf = 9.4 cm/s). However, at 



small /, P(I) develops an upward trend which becomes 
increasingly more evident as the flow velocity decreases. 

In this letter, we report the results of event-driven 
simulations of a system of inelastic, monodispersed hard 
disks in the experimental geometry of Ref. 11]. Our sim- 
ulations provide clear evidence of an increasing propor- 
tion of collisions with small impulses as the flow velocity 
is decreased. This results in the formation of a plateau 
in P(I) as the minimum flow velocity for sustained flow 
is approached. In addition, we observe the formation of 
clusters of disks which collide "frequently" and are rem- 
iniscent of the "collaose strings" observed in freely cool- 
ing granular matter[l2j. We present a model calculation 
which strongly suggests that the increase of small impulse 
events is associated with the growth of these clusters. 

Simulations The grain dynamics used in the simula- 
tions are as in Ref. [13| (momentum is conserved and at 
each interparticle collision the energy loss is proportional 
to (1 — fi 2 ) where /x is the coefficient of restitution). To 
ensure that the pressure is independent of the height the 
side walls must absorb some vertical momentum, there- 
fore we impose the condition that collisions with the walls 
are inelastic in the tangential direction. The flow velocity 
is controlled similarly to the experiments, by adjusting 
the width of the hopper opening. However, as we wish 
to observe the system over many events, particles exiting 
the system at the bottom must be replaced at the top to 
create sustained flow. This necessitates the introduction 
of a probability of reflection p at the bottom of the hop- 
per (this would be equivalent to the presence of a sieve in 
experiment). In our case, it provides another parameter 
with which we can tune the flow velocity. Typically, our 
simulations were done on systems of 500 disks, with /i = 
0.9 and p = 0.4. The simulation was run for 2 x 10 6 events 
for each flow velocity, with 1.5 x 10 5 discarded initially 
to allow the system time to reach steady state. 

Simulation Results The physical quantity that is 
most closely related to P(f) in our flowing hard-disk sys- 
tem is the distribution of impulses transferred at each 
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FIG. 1: Impulse distribution scaled to the average im- 
pulse, <I>, with flow velocity decreasing from the bot- 
tom curve to the top curve. The inset shows the ratio 
of P(Ipeak) to P(I m in) as a function of the flow velocity. 



collision, P(I). Unlike the experiments done in Ref. 
where the measurements were made only at the wall, 
we measure the magnitude of the momentum transfer 
at each collision event (note that the impulse and flow 
velocity are measured in units such that disk diameter 
d, disk mass m and acceleration due to gravity g are all 
equal to one). We observe an exponential form of P(I) 
at impulses larger than the average impulse (Fig. ^) for 
the full range of flow rates. After scaling the impulse 
by the average impulse, < I >, the curves collapse onto 
each other for large impulses. The small impulse be- 
havior of the distribution, however, changes markedly 
with changes in the flow velocity. As the flow velocity 
decreases, the height of the distribution at a minimum 
impulse (I m in — 0.0075 < I > for this distribution) in- 
creases and P(I) begins to develop a plateau at small /. 
If the ratio between P(I) at the peak impulse to P(I) 
at the minimum impulse l m in is calculated and plotted 
against the flow velocity, a dramatic increase is observed 
(see inset to Fig. Extrapolating to small flow veloc- 
ities indicates that this ratio becomes one at a non-zero 
flow rate. 

It should also be noted that the height of the distribu- 
tion at / = makes a sharp upturn which becomes more 
pronounced as the flow velocity decreases (the upturn is 
not present for the fastest flow velocity). While we do 
not fully understand the reasons for this phenomenon, 
we believe it is not the same upward trend visible in 
the experiments. That feature was a smooth continua- 
tion of the large impulse curve, while in our distribution 
a peak is seen. In addition, our upturn occurs only at 



1 = 0, regardless of the bin width used in constructing 
the distribution. Such small impulses are not likely to 
be distinguishable in experiment due to the resolution 
limit imposed by the transducer. For the remainder of 
our analysis, we will consider the impulse distribution 
without the 1 = point. 

The basic form of the impulse distribution, a peak at 
a finite value of / and an exponential tail can be under- 
stood on the basis of a flow of uncorrelated hard disks. 
If the disks were uncorrelated, the impulse distribution 
would be a convolution of the individual momentum (ve- 
locity) distributions. Since there is an average flow ve- 
locity, this would give rise to a peak in the distribution. 
The exponential tail would be observed if the individual 
velocity distributions also had exponential tails. The ve- 
locity distribution observed in our simulations can, to a 
first approximation, be described by Gaussians with ex- 
ponential tails. This argument would then suggest that 
the exponential tail arises from uncorrelated particles and 
is a consequence of the shape of the velocity distribution. 
Viewed from this perspective, the lack of change in shape 
of P(I) at large / as the flow approaches jamming is not 
surprising. We will discuss this further in the context of 
a simple model to be presented below. 

The change in P(I) at small I is reminiscent of the uni- 
versal jamming scenario |^ although the direction of the 
change (filling up at small impulses) is different from the 
one observed in Lennard Jones systems and foams where 
the probability of finding small forces decreases as the 
system nears the jamming transition. Nevertheless, the 
idea that changes in the distribution at small impulses 
or forces could be related to the appearance of structures 
akin to force chains is intriguing. 

To explore a possible connection between the changes 
in P(I) and the appearance of spatial inhomogeneities, 
we considered a question first asked in studies of inelastic 
collapse in freely cooling granular gases [l^. "How many 
collisions does a given grain undergo in a fixed number of 
events?" We define a minimum frequency of collision for 
a given particle as u>o = ^g= (where r is the average time 
between events for a given flow velocity) and then iden- 
tify those particles undergoing collisions with a frequency 
> ujq as frequently colliding. This allows us to define a 
minimum number of collisions a frequently colliding par- 
ticle must undergo in a given number of events. We then 
construct an image of our system at every 1000 events, 
and color all disks satisfying the criteria (see Fig. [5Ji,b). 
As we decrease the flow velocity, the frequently colliding 
particles form increasingly larger linear clusters (compare 
Fig. UK) where Vf = 2.03 in our units or 35.6 cm/s and 
Fig. 03, where Vf = 0.96 or 16.9 cm/s). These ID struc- 
tures observed in our simulation are reminiscent of the 
transient solid chains postulated by the hydrodynamic 
model of Ref.Q. Comparisons of the impulse distribu- 
tion of the frequently colliding particles and the impulse 
distribution of the remaining "rarely colliding" particles 
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FIG. 2: (a) (b) Sample image of simulation for Vf — 
2.03 (a) and for Vf — 0.96 (b). (c) (d) Impulse distribu- 
tions for all particles, frequently colliding particles and 
rarely colliding particles for vf = 2.03 (c) and for vt = 
0.96 (d). Note that these distributions are unnormalized. 



(Fig. |21 c,d) reveals that it is the contribution from the 
frequently colliding particles which causes the height of 
the total impulse distribution at I m i n to increase in the 
manner seen in Fig. The large impulse behavior ap- 
pears to be dominated by the rarely colliding particles. 
These observations are relatively insensitive to changes 
in loq provided it is larger than a minimum value. 

Impulse Distribution Model A possible connection be- 
tween the shape of the impulse distribution and the de- 
velopment of linear clusters of frequently colliding parti- 
cles can be drawn by investigating the following model. 
Consider a one dimensonal cluster of particles which are 
all moving with the same velocity. Now if another par- 
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FIG. 3: (a) Particle of speed v incident on ID cluster of 
S stationary particles, (b) Resulting form of P(I,v, S). 



tide, travelling at some speed v relative to the cluster 
(see Fig. collides inelastically with one end of the 
chain, the impulse associated with that collision will be 
(1 — e)v where e is related to the coefficient of restitution 
by e = [13 • If the chain was only comprised of one 
grain, the impulse distribution P(I, v, 1) (given an initial 
incoming speed v) would be a single spike of height 1 lo- 
cated at (1 — e)v. For two particles, P(I,v,2) would be 
a bar of height 1 with limits at (1 — e) 2 v and (1 — e)v. 
Continuing this argument for clusters containing S par- 
ticles, then as S becomes very large, the leftmost limit of 
P(I, v, S) approaches zero (Fig. Oh). 

Given a distribution of speeds P(v) for the incident 
particle, and a distribution of cluster sizes P(S), the total 
impulse distribution P(I) is: 



P(I) = J dSdvP(v)P(S)P{I,v,S) 



(1) 



If the cluster size distribution falls off sharply, then 
the shape of P(I) will essentially reflect the shape of 
P(v). If, however, the cluster size distribution becomes 
broad, then the small impulse end of P(I) will flatten 
out reflecting the nature of P(I, v, S) for large S. 

As detailed previously, the criterion that we use to 
identify the clusters in our simulations is that of "fre- 
quent" collisions. Due to the inelastic nature of the col- 
lisions, the velocities of these particles become highly 
correlated^ and they can be idealized as clusters of 
disks moving with the same velocity. This velocity corre- 
lation within spatial clusters of particles has been directly 
observed in granular surface flows 8]. Since we observe 
these clusters to be linear and growing in size with de- 
creasing flow rates, it is plausible that the changes in P(I) 
observed in our simulations is related to a change in dis- 
tribution of the size of the clusters of frequently-colliding 
disks. In order to put this conjecture on a firmer footing, 
we calculated the impulse distribution from our model 
using forms of P(v) and P(S) which provide good fits to 
our simulation data. The form of P(v) that we use is a 
Gaussian with an exponential tail. This form is represen- 
tative of the simulation results for P(v) with the center 
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suggests that the observed changes in P(I) are associ- 
ated with the growth of clusters of frequently colliding 
particles and as S max diverges, the ratio of P{I m in) to 
P(Ipeak) approaches unity. The tail of the P(I) distribu- 
tion is controlled by the shape of P(v) as we have verified 
within the model. 

The picture which seems to be emerging from our sim- 
ulations is that of increasingly larger scale spatial hetero- 
geneities developing as the system approaches jamming. 
Since the clusters are essentially the same as the ones 
identified in freely cooling granular matter, their origin 
lies in the dissipative nature of the medium. The het- 
erogeneities reflect strong velocity correlations of grains 
and leave a distinctive signature in the impulse distribu- 
tion. Whether the clusters that we have identified are in- 
deed incipient force chains remains to be verified. If this 
connection can be established, for example through the 
calculation of stress correlations in the flowing medium, 
then our analysis will provide a natural connection be- 
tween P(f) and force chains. Moreover, it would a in- 
dicate that jamming occurs through the formation of 
system spanning clusters of grains whose velocities are 
strongly correlated. 

Acknowledgements We thank N. Menon, S. Tewari, 
N. Easwar and S.R. Nagel for many helpful discussions. 
AF, BF and BC acknowledge support from NSF through 
grant No. DMR 0207106, and AF acknowledges sup- 
port from the Natural Sciences and Engineering Research 
Council, Canada. 



1.0 0.5 1.0 1.5 
!/<!> 



2.0 



FIG. 4: (a) Simulation results for SP(S) for vary- 
ing flow rates. (b) Results of the proposed model 
for P{I) with S m ax increasing from the bottom 
curve to the top curve. The inset shows the ra- 
tio Of P{Ipeak) tO P(Irnin) aS & function of 1/ S max . 



of the Gaussian reflecting the average flow velocity. The 
observed cluster size distribution is consistent with the 
form P(S) = cxp(— S/S max )/S with S max increasing as 
the flow velocity decreases. Fig.QJi shows P(S) obtained 
from our simulations for three different flow rates. Note 
that a change in loq will produce a change in S max for a 
given flow velocity but will not alter the shape of P(S) 
nor the trend of increasing S ma x with decreasing flow 
velocity. 

Using these forms of P(v) and P(S), we find that the 
P(I) obtained from the simple model has a form which 
is remarkably similar to the measured P(I), as shown in 
Fig-EJ). The control parameter in the model is S max and 
the contribution at small impulses increases with S max 
as is evident from the inset to Fig. The correspon- 
dence between the model and our simulations strongly 
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